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Abstract  Marine  phytoplankton  and  associated  organic  materials  absorb  a  substantial  quantity  of  solar 
shortwave  energy  penetrating  the  upper  ocean.  Most  of  this  absorbed  energy  is  lost  as  heat  and  thereby 
contributes  to  the  warming  of  near-surface  waters.  Here  we  examine  this  biothermal  feedback  effect  on  upper 
ocean  physics  and  air-sea  energy  exchange  using  a  fully  integrated  ocean-atmosphere-biological  modeling 
system.  Our  model  simulations  show  that  a  local  phytoplankton  bloom  may  impact  upper  ocean  physics  in  such 
a  way  as  to  promote  the  spatiotemporal  persistence  of  the  bloom  itself  within  a  semi-enclosed  coastal 
embayment.  This  is  accomplished  primarily  via  enhanced  thermal  stratification  that  promotes  vertical  stability 
and  more  efficient  utilization  of  macronutrients.  Modulations  of  wind  stress  patterns  due  to  perturbations  in  the 
local  surface  pressure  gradients  also  arise  as  a  result  of  the  simulated  biothermal  warming  of  surface  waters.  The 
model  evidence  suggests  that  the  observed  persistence  of  phytoplankton  blooms  in  the  northern  Monterey 
Bay,  California,  may  be  enhanced  by  similar  synergistic  interactions  between  ocean  biology  and  physics. 
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1.  Introduction 

A  fundamental  concept  in  biological  oceanography  is  that  the  physical  processes  of  the  oceans  largely 
determine  the  spatiotemporal  variability  of  phytoplankton  abundance  and  productivity.  For  example,  the 
wind-driven  coastal  divergence  of  surface  waters  along  the  eastern  ocean  margins  results  in  the  upwelling  of 
nutrient-rich  deeper  waters  that  stimulate  microalgal  growth  [Chavez  and  Messie,  2009;  Walsh,  1988],  The 
main  features  of  global  phytoplankton  abundance  evident  in  synoptic  satellite  data  may  indeed  be  broadly 
explained  in  this  context  of  geophysical  forcing  [see  Longhurst,  1998], 

However,  the  abundance  of  marine  phytoplankton  belies  the  inefficiency  of  photosynthesis  as  a 
photochemical  process:  most  of  the  light  energy  absorbed  by  algal  pigments  is  lost  to  the  surrounding 
environment  as  heat  [Morel,  1978;  Bannister  and  Weidemann,  1984;  Morel,  1988].  Accordingly,  it  has  been 
theorized  that  marine  phytoplankton  have  the  potential  capacity  to  modulate  the  heating  of  the  upper  ocean 
due  to  the  optical  properties  of  microalgal  pigments  [Lewis  et  at.,  1983;  Morel,  1988].  Furthermore,  the  high 
turnover  rate  of  algal  biomass  (~  2  days  [Falkowski,  1994])  generates  more  temporally  persistent  nonliving 
organic  matter  that  absorbs  substantial  quantities  of  solar  energy  in  the  surface  ocean  [Bricaud  et  a!.,  1981]. 
Thus,  the  aggregate  biothermal  impact  of  phytoplankton  and  associated  organic  constituents  upon  the  heat 
budget  of  the  upper  ocean  is  not  negligible  [Morel  and  Antoine,  1 994]. 

Quantitatively  resolving  the  biothermal  feedback  effect  upon  oceanographic  processes  is  nonetheless  difficult. 
An  observational  study  would  ideally  require  a  control  so  that  identical  atmospheric  and  physical 
oceanographic  conditions  may  be  experienced  with  and  without  the  additional  optical  attenuation  provided  by 
phytoplankton  and  associated  organic  constituents.  Given  that  this  approach  is  not  feasible,  most  of  the 
progress  toward  understanding  the  biothermal  impacts  has  been  accomplished  via  ocean  models  where  such 
an  experiment  may  be  performed  within  the  simulations  [Anderson  et  at.,  2007;  Cahill  et  at.,  2008;  Manizza  et  al., 
2008;  Oschlies,  2004;  Rochford  et  al.,  2001 ;  Sweeney  et  al.,  2003;  Wu  et  al.,  2007].  More  recently,  two-way  coupled 
ocean-atmosphere  numerical  models  suggest  that  this  biothermal  effect  is  significant  on  global  and  climatic 
scales  [Patara  et  al.,  2012]  as  well  as  locally  and  on  much  shorter  timescales  [Jolliff  et  al.,  2012a].  Two-way 
coupled  ocean-atmosphere  modeling  systems  better  approximate  the  air-sea  exchange  of  thermal  energy  as 
constrained  by  conservation;  however,  continuing  dynamical  impacts  back  upon  the  surface  biota  cannot  be 
assessed  without  an  additional  biological  model.  In  this  paper,  this  deficiency  is  addressed  by  integrating  a 
biological  model  into  a  numerical  ocean-atmosphere  modeling  system.  This  allows  for  not  merely  an 
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(A)  COAMPS-B  [SI] 


assessment  of  how  biothermal  effects  may  impact  the 
physics  but  also  a  further  analysis  of  how  these  physical 
impacts  may,  in  turn,  impact  the  biology. 


(B)  COAMPS-B  [S2] 


2.  Methods 

The  Naval  Research  Laboratory  (NRL)  has  developed  a 
fully  integrated  ocean-atmosphere-biological  numerical 
modeling  system  based  on  the  Coupled  Ocean- 
Atmosphere  Mesoscale  Prediction  System  (COAMPS®). 
COAMPS  is  a  nested,  relocatable,  and  two-way  coupled 
ocean-atmosphere  modeling  system  that  is  presently 
used  for  ocean-atmosphere  forecasting  by  the  U.S.  Navy. 
The  nonhydrostatic  atmospheric  COAMPS  model 
component  [ Hodur ,  1997]  is  the  Navy's  operational 
mesoscale  meterological  forecasting  system.  The 
hydrostatic  Navy  Coastal  Ocean  Model  (NCOM)  [ Barron 
et  at.,  2004]  serves  as  the  oceanic  component  of 
COAMPS.  This  modeling  component  was  executed  in 
data-assimilative  mode  via  the  Navy  Coupled  Ocean 
Data  Assimilation  system  [Cummings,  2005]  for  a 
1  month  spin-up  (beginning  on  1  April  2008)  and  then  in 
non-data-assimilative  (or  "free-run"  mode)  for  the 
remaining  3  months  of  the  simulation  (ending  31  July 
2008).  The  atmospheric  and  oceanic  model  coupling  was 
designated  via  the  uppermost  oceanic  model  grid  cell 
temperature  and  the  lowest  grid  cell  atmospheric  model 
variables  (temperature,  humidity,  wind  velocity, 
pressure,  and  radiative  fluxes).  At  a  6-min  coupling 
interval,  bulk  fluxes  of  thermal  energy  exchange  were 
calculated  following  the  Coupled  Ocean-Atmosphere 
Response,  version  3,  scheme  [ Fairali  et  at.,  1996]. 

The  integration  of  a  four-component,  nitrogen-based 
biological  model  into  COAMPS  (Figure  1)  occurs  first  via 
the  conventional  means:  biological  state  variables 
(phytoplankton  biomass,  nitrate  nitrogen,  ammonium 
nitrogen,  and  particulate  detritus  nitrogen)  are  physically 
forced  by  the  advection  and  diffusion  resolved  by  the 
ocean  circulation  model  (NCOM).  The  physical  ocean 
model,  however,  must  also  be  informed  about  how  to 
attenuate  solar  energy  to  meet  its  requirement  for  solar 
heating  rate  computations.  The  ocean  biology  model 
provides  this  information.  Biological  state  variables  are 
used  to  estimate  the  bio-optical  attenuation,  and  these 
terms  are  provided  back  to  the  physical  ocean  model.  Thus,  changes  in  the  biological  state  variables  have  the 
potential  to  impact  the  simulated  thermodynamic  upper  ocean  processes.  The  specific  details  of  the 
biological  model  and  its  interface  with  the  physical  models  are  provided  in  Appendix  A. 


Figure  1.  Schematic  of  the  COAMPS-B  modeling  sys¬ 
tem.  Arrows  and  text  summarize  the  communication 
between  the  modeling  components,  (a)  The  model 
schematic  for  COAMPS-B  with  the  biothermal  feedback 
disabled  (SI).,  (b)  Biothermal  feedback  as  the  transfer  of 
vertical  optical  attenuation  coefficients  the  physical 
model  uses  to  determine  the  penetration  of  solar 
shortwave  radiation  into  the  surface  ocean.  This  model 
(Figure  1b)  is  used  for  simulation  2  (S2). 


To  isolate  the  potential  impact  of  the  biothermal  feedback  in  the  model,  two  simulations  are  performed:  (SI)  the 
feedback  from  the  biology  is  disabled  (Figure  la);  and  (S2)  the  feedback  from  the  biology  to  the  physics  is  active 
(as  in  Figure  1  b).  In  the  former  simulation  (SI),  the  physical  ocean  model  (NCOM)  presumes  a  constant  set  of 
optical  attenuation  terms  based  on  the  work  of  Paulson  and  Simpson  [1977].  These  computations  are  an 
approximation  of  the  Jerlov  water  types  [Jerlov,  1 976],  and  they  are  often  used  in  physical  ocean  models  [Kara 
etai.,  2004;  Rochford  et  at.,  2001].  Whereas  these  terms  may  constitute  a  reasonable  approach  in  some  cases,  the 
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main  point  of  criticism  is  that  they  are  not 
dynamic,  i.e.,  the  optical  attenuation  terms  are 
invariant  and  so  they  cannot  mimic  the  dynamic 
spatiotemporal  changes  in  ocean  bio-optical 
properties  that  may  occur  in  many 
coastal  settings. 


3.  Study  Site 

The  COAMPS-B  (COAMPS  with  a  Biological 
module)  nested  domain  is  established  over 
Monterey  Bay,  California  (herein  "MB";  Figure  2). 
It  is  important  to  note  that  verification  and 
validation  of  the  COAMPS  forecasting  system 
may  be  found  elsewhere  [Doyle  et  al.,  2009; 
Eliasen  et  al.,  2005;  Small  et  al.,  201 2].  Here  the 
focus  is  placed  on  the  modeling  system's 
sensitivity  to  changes  in  the  oceanic  shortwave 
attenuation  that  arise  as  a  result  of  the  dynamic 
surface  ocean  phytoplankton  variability  that  is 
simulated  by  an  ocean  biological  model. 

From  March  through  November,  the  coastal 
region  surrounding  MB  is  subject  to  sustained 
northwesterly  winds  resulting  in  coastal 
upwelling.  Periodic  slackening  of  these 
upwelling  winds  are  referenced  as  "relaxation" 
events  [ Breaker  and  Broenkow,  1 994;  Pennington 
and  Chavez,  2000].  During  upwelling  winds,  cold 
water  filaments  north  and  south  of  the  bay  arise 
from  the  coastal  divergence  and  a  strong 
equatorward  flow  of  surface  waters  is  often 
observed  across  the  mouth  of  the  bay. 
Divergence  of  the  southward  flow  toward  the 
southern  portion  of  the  bay  brings  about  a  semicyclonic  circulation  within  the  bay  [see  Ramp  et  al.,  2005].  The 
intense  surface  currents  (a  surface  "jet")  across  the  mouth  of  the  bay  prevent  egress  of  the  surface  waters  out  its 
northwestern  periphery.  The  Santa  Cruz  Mountains  immediately  to  the  north  of  this  area  effectively  block  the 
prevailing  northwesterly  winds.  This  combination  of  physical  conditions  sets  up  a  relatively  calm  and  retentive 
regime  for  the  northeastern  waters  of  Monterey  Bay.  Hence,  this  area  is  referred  to  as  the  "upwelling  shadow" 
[Graham  and  Largier,  1 997].  The  retentive  nature  of  the  surface  flow  subjects  these  waters  to  elevated  surface 
warming  and  enhanced  biological  productivity  such  that  satellite  detection  of  sea  surface  temperature  (SST) 
and  surface  chlorophyll  routinely  depict  elevated  signals  in  northern  Monterey  Bay  [e.g.,  Ryan  et  al.,  2010]. 

Herein  COAMPS-B  is  used  to  simulate  the  physical  and  biological  processes  that  occur  in  Monterey  Bay's 
upwelling  shadow.  The  objective  is  to  isolate  the  potential  for  biothermal  effects  to  impact  the  physical  and 
biological  processes  that  occur  within  this  semi-enclosed  embayment.  In  section  4.1,  the  differences  in  the 
physical  simulations  are  examined  when  biothermal  effects  are  accounted  for  in  the  modeling  system.  In 
section  4.2,  the  potential  for  synergy  between  the  interacting  biological  and  physical  fields  are  examined,  and 
section  4.3  examines  in  situ  data  in  light  of  these  analyses.  The  in  situ  data  were  collected  in  Monterey  Bay 
during  June  2008  as  part  of  an  NRL  field  program;  additional  details  may  be  found  in  Jolliffet  al.  [2012b]. 

4.  Results 

4.1.  Biothermal  Impact  on  the  Simulated  Physics 

The  differences  between  the  simulations  demonstrate  the  potential  biothermal  perturbation  provided  by  the 
presence  of  phytoplankton  in  the  surface  ocean.  For  example,  a  longitudinal  cross  section  through  the 


Longitude  (°W) 


Figure  2.  (a)  The  COAMPS  atmospheric  model  nests  beginning 
at  27  km  horizontal  resolution  and  scaling  down  to  the  inner 
nest  of  3  km  resolution  over  the  Monterey  Bay,  California,  area, 
(b)  The  innermost  atmospheric  nest  (1  km)  overlaid  the  inner 
ocean  model  nest  (500  m  horizontal  resolution). 
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COAMPS-B  model  results  for  the  ocean 
and  atmosphere  (Figure  3)  reveals  an 
upper  surface  ocean  (<  10  m)  warming 
and  enhanced  thermal  stratification  in 
S2  that  is  not  present  in  SI ,  particularly  in 
the  northern  portion  of  the  transect 
(north  of  36.8°N).  The  simulations  are 
two-way  coupled  between  the 
atmospheric  and  oceanic  model 
components.  The  biothermal 
perturbation  present  in  the  ocean  is 
thereby  transmitted  to  the  simulated 
lower  atmosphere.  The  lower  atmosphere 
in  the  northern  portion  of  the  bay  is  as 
much  as  ~1°C  warmer  in  S2  (Figure  3). 

These  coupled  ocean-atmosphere 
simulation  results  are  very  similar  to 
what  was  reported  in  Jolliffetal.  [2012a]. 
In  the  earlier  work,  satellite  data  rather 
than  an  integrated  biological  model 
determined  the  oceanic  shortwave 
penetration  depth.  These  simulations 
also  suggest  that  the  additional  thermal 
energy  retained  by  phytoplankton  near 
the  surface  results  in  elevated  turbulent 
heat  flux  transfers  from  the  ocean  to  the 
atmosphere.  These  differences  are  acute 
during  the  relaxation  periods  of  the 
prevailing  northwesterly  winds 
(Figure  4a),  which  were  mimicked  by  the 
COAMPS  atmospheric  model 
component  (Figure  4b)  from 
approximately  13  to  15  June. 


Figure  3.  (a)  A  longitudinal  cross  section  of  the  simulated  (top)  atmospheric 

temperatures  and  (bottom)  ocean  temperatures  for  simulation  SI.  (b)  The  Differences  in  the  mean  SST  fields 
same  results  for  S2.  Both  results  are  from  14  June  2008  at  1200  GMT.  immediately  before  (10-12  June)  and 

during  (13-15  June)  the  midmonth 

relaxation  event  reveal  a  substantial  difference  in  surface  warming  (Figures  5a  and  5b).  This  heating  difference 
may  also  be  shown  by  computing  the  spatially  averaged  SST  in  northern  MB  (north  of  36.8°N  and  west  of 
—  122.1°W,  as  shown  in  Figures  5a  and  5b)  and  calculating  the  difference  between  the  respective  simulation 
means  at  the  hourly  simulation  output  increment  (S2  —  SI;  Figure  5c).  By  300  h  into  the  June  simulation,  the 
mean  SST  differences  begin  to  peak  at  ~1°C.  These  temperature  differences  are  transmitted  to  the  lower 
atmosphere  overlying  the  northern  bay  as  a  temperature  difference  also  approaching  ~1°C  (Figure  6a).  This 
warmer  air  has  a  greater  capacity  to  retain  water  vapor,  and  the  simulated  relative  humidity  is  thereby  reduced 
downward  by  as  much  as  8%  (Figure  6b).  This  relative  humidity  decrease  obscures  the  overall  modest  increase 
in  total  water  vapor  content  for  S2  (data  not  shown).  This  increased  heat  and  moisture  content  in  the  lower 
atmosphere  of  S2  is  principally  due  to  the  increase  in  latent  heat  flux  resolved  by  the  simulation.  The  spatially 
averaged  latent  heat  flux  differences  peak  as  high  as  1 5  W  mT2  after  300  h  (Figure  6c). 


The  simulated  presence  of  phytoplankton  warms  the  upper  ocean,  increases  the  upper  ocean  thermal 
stratification,  and  increases  the  temperature  of  the  overlying  lower  atmosphere.  These  simulated  effects  are 
conspicuous  during  the  wind  relaxation  periods  when  local  effects  prevail,  as  opposed  to  the  dominant 
forcing  of  the  regional-scale  wind  stress  patterns.  From  an  energy  conservation  point  of  view,  the  heat 
retentive  capacity  of  the  phytoplankton  in  the  upper  ocean  effectively  traps  a  greater  amount  of  solar 
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(A) 
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Figure  4.  (a)  Measured  wind  velocities  (ms  1 )  at  Monterey  Bay  Aquarium  Research  Institute  (MBARI)  mooring  Ml  from  1  to 
30  June  2008;  (b)  the  corresponding  wind  velocities  from  COAMPS. 


(A)  SST  Bias  10-12  June  2008 


(B)  SST  Bias  13-15  June  2008 
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Figure  5.  (a)  The  difference  in  the  temporal  means  in  SST  (°C)  between  S2  and  SI  for  the  period  10-12  June,  (b)  The 
difference  in  the  temporal  means  in  SST  between  S2  and  SI  for  the  period  13-15  June,  (c)  The  hourly  difference  in  the 
mean  SST  fields  (S2  —  SI)  for  June.  The  means  are  drawn  from  the  area  of  northern  MB  indicated  in  Figures  5a  and  5b. 
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Hours  from  1  June  00  GMT 


Figure  6.  (a)  The  difference  in  the  temporal  means  in  2  m  air  temperature  (°C)  between  S2  and  SI  for  the  period  13-15 
June,  (b)  The  difference  in  the  temporal  means  in  the  relative  humidity  (%)  between  S2  and  SI  for  the  period  13-15 
June,  (c)  The  hourly  difference  in  the  mean  latent  heat  flux  fields  (S2  -  SI;  W  m-2)  for  June.  The  means  are  drawn  from  the 
area  of  northern  MB  indicated  in  Figures  5a  and  5b. 


shortwave  near  the  air-sea  interface.  This  energy  increase  is  readily  fluxed  back  into  the  lower  atmosphere. 
There  are  two  additional  simulated  consequences  to  the  biothermal  perturbation  resolved  here:  (1)  the 
changes  in  fluid  densities  and  energy  fluxes  appear  to  have  some  influence  on  simulated  local  wind  patterns 
and  upper  ocean  circulation  during  the  wind  relaxation  period;  and  (2)  the  increased  near-surface  ocean 
stratification  appears  to  amplify  the  surface  biological  productivity. 

The  first  of  these  effects  is  minor  but  cannot  be  summarily  dismissed  as  negligible.  For  example,  the  lower 
level  heating  and  moisture  increase  in  S2  introduces  a  -0.04  hPa  lowering  of  sea  level  pressure  over  the 
northern  bay  (Figure  7a).  Nonetheless,  this  modest  decrease  is  enough  of  a  local  deviation  in  the  zonal 
pressure  gradient  to  comparatively  accelerate  the  mean  westerly  component  of  the  wind  velocity 
by -0.3  ms-1  (Figure  7b).  The  simulated  instantaneous  wind  velocity  differences  approach  1.0  ms-1  in  the 
center  of  Monterey  Bay  on  14  June  2008. 

The  magnitude  of  the  wind  velocity  differences  is  not  large  since  the  overall  wind  forcing  is  diminished 
during  the  relaxation  period.  The  salient  point  is  that  during  periods  when  local  conditions  within  Monterey 
Bay  dominate  the  physical  dynamics,  as  when  the  regional  wind  patterns  driven  by  the  larger-scale  pressure 
gradients  are  diminished,  the  air  and  sea  planetary  boundary  layer  circulation  resolved  by  the  two 
simulations  become  incoherent.  For  example,  the  SI  and  S2  north/south  surface  water  transports  through 
the  bay  (through  36.8°N)  are  highly  correlated  until  the  13  June  relaxation  period  when  the  transports 
become  significantly  out  of  phase  (Figure  8,  solid  line).  This  episode  corresponds  to  a  mismatch  in  the 
simulated  westerly  wind  stress  magnitudes  over  the  bay  (Figure  8,  thick  line).  The  eddy  kinetic  energy  of  the 
surface  currents  (computed  as  [u2  +  v1V2)  is  also  diminished  in  S2  during  the  relaxation  period  (Figure  8, 
dashed  line). 
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Figure  7.  (a)  The  difference  in  the  temporal  means  for  sea  level  pressure  (hPa)  between  S2  and  SI  for  the  period  13-15 
June,  (b)  The  difference  in  the  temporal  means  of  the  westerly  component  of  the  winds  (m  s_1)  between  S2  and  SI  for 
the  period  13-15  June. 


4.2.  Synergy  Between  Biological  and  Physical  Fields 

Simulated  upper  ocean  chlorophyll  distributions  are  spatially  coherent  to  those  of  SSTin  the  northern  portion 
of  MB  during  the  wind  relaxation  regime  (Figure  9).  The  amplification  of  both  chlorophyll  and  SST  from  SI  to 
S2  suggests  a  synergy  between  the  two  variables.  The  most  obvious  mechanism  is  that  the  enhanced  thermal 
stratification  in  S2  promotes  enhanced  phytoplankton  growth  near  the  surface,  which  in  turn,  promotes  a 
further  enhancement  of  thermal  stratification  via  the  biothermal  feedback  effect.  The  simulations  quantify 
this  potential  synergy  in  northern  MB  as  approximately  0.8°C  in  increased  SST  and  2.9  mmol  C  rrT3  additional 
surface  phytoplankton  biomass  on  14  June  (1200  GMT).  These  are  the  differences  in  spatial  simulation 
averages  from  northern  MB;  differences  in  average  quantities  are  referred  to  herein  as  a  "bias."  The  biomass 
bias  is  particularly  significant  when  examined  in  light  of  the  SI  average  biomass  quantity,  i.e.,  the  ratio  of  the 
biomass  bias  [mean(SI)  —  mean(S2)]  to  the  initial  simulation  mean  [mean(SI)].  The  biothermal  effect 
increases  the  simulated  surface  biomass  by  27%. 


Figure  8.  The  linear  correlation  between  north/south  surface  water 
transport  (m3s_1)  to  10m  depth  along  36.8°N  from  -122.1°W  to  the 
coastal  boundary  is  shown  by  the  thin  line.  The  thick  line  is  the  ratio  of 
the  S2  to  SI  wind  stress  (zu(S2)/  tw(S1  )  —  1 .0);  the  dashed  line  is  the  ratio 
of  surface  current  eddy  kinetic  energy  EKE  (EKE(S2)/EKE(S1)  -  1.0), 
where  EKE  =  ([u2  +  v2]/2).  EKE  and  zu  values  were  first  smoothed  with  a 
running  24  h  (24  point)  average. 


Thermal  upper  ocean  stratification  within 
northern  MB,  quantified  as  AT  from  the 
surface  to  22  m  depth  (A r22),  is  biased 
upward  by  ~1°C  in  S2  during  the  simulated 
wind  relaxation  period  (Figure  10).  SST  and 
surface  chlorophyll  are  similarly  elevated 
(by  an  average  of  0.8°C  and  0.7  mg  chi  rrT3, 
respectively;  Figure  10).  The  simulations  do 
not  feature  variable  carbon-to-chlorophyll 
ratios  (see  section  A1);  these  differences  in 
chlorophyll  are  directly  forced  by  differences 
in  simulated  phytoplankton  biomass. 

Indeed,  the  relative  simplicity  of  the 
biological  simulation  (four  nitrogen-based 
compartments)  enables  a  conceptually 
simple  contrast  between  the  simulated 
nitrogen  budgets.  Defining  the  organic 
nitrogen  as  that  within  the  phytoplankton 
and  detritus  reservoirs,  and  the  inorganic 
nitrogen  as  the  nitrate  and  ammonium 
reserves,  the  total  amount  of  organic 
nitrogen  within  the  upper  1 0  m  of 
northern  MB  increases  by  an  average  of 
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Figure  9.  (a)  SI  surface  chlorophyll  (mg  chi  m  3)  on  14  June 
SI  on  14  June  1200  GMT;  (d)  the  same  for  S2. 
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GMT;  (b)  the  same  for  simulation  S2.  (c)  SST  for  simulation 


36.6N 
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(D)  S2  SST,  14  June  1200  GMT 
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1 0%  in  simulation  S2  with  a  peak  increase  of  23%  during  the  wind  relaxation  period  (Figure  1 1 ).  By  the  end  of 
June,  S2  has  approximately  10%  less  inorganic  N  in  the  surface  layer  than  SI  (Figure  11).  The  enhanced 
thermal  stratification  in  S2  has  a  tendency  to  inhibit  turbulent  vertical  diffusion,  which  would  otherwise  tend 
to  homogenize  the  vertical  distribution  of  scalar  quantities.  The  additional  biological  benefit  is  that 

phytoplankton  are  more  restricted  to  the 
near-surface  layer — sustaining  continuing 
exposure  to  available  light  and  improving 
the  utilization  of  inorganic  nitrogen  for 
growth.  This  is  particularly  advantageous 
when  light  is  more  limiting  to  growth  than 
available  macronutrients. 

Simultaneously,  the  increased  thermal 
stratification  in  S2  restricts  the  egress  of 
thermal  energy  out  of  the  near-surface 
layer  (<  1 0  m  depth),  and  this  also 
impacts  the  thermal  energy  budgets  in 
addition  to  the  nitrogen  budgets.  The 
respective  diffusive  differences  tend  to 
accelerate  the  rate  of  warming  in  S2 

Figure  10.  The  hourly  difference  in  mean  SST  (°C)  for  northern  MB  is  during  periods  wherein  thermal 
depicted  by  the  thick  line.  The  difference  in  A T22  (°C;  an  index  of  stra-  stratification  prevaNs  in  both  simulations, 
tification,  as  explained  in  the  text)  is  shown  by  the  dashed  line,  and  the  . 

chlorophyll  difference  (mg  m-3)  is  shown  by  the  gray  line.  In  all  cases,  For  sample,  the  S2  depth-averaged 

the  difference  is  expressed  as  S2  —  SI ;  positive  differences  indicate  temperature  in  the  upper  10.7  m  from  1 2 

greater  heat  and  biomass  in  the  surface  layers  of  simulation  S2.  June  (2000  GMT)  to  1 5  June  (0600)  in 
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Differences  in  Total  Nitrogen  Distributions:  (S2  -  S1)/S1 


Figure  11.  The  thin  line  is  the  difference  (as  a  %)  between  organic 
nitrogen  (as  defined  in  the  text)  in  S2  versus  SI  in  the  upper  10  m  for 
northern  MB.  The  total  quantity  of  organic  nitrogen  in  each  simulation 
is  used  to  calculate  the  difference  as  (S2  -  SI )/S1  x  100.  The  thick  line  is 
the  inorganic  nitrogen  difference  (as  a  %). 


northern  MB  increased  at  a  rate  of  0.24° 
Cd~\  Given  the  expression: 


AT  _  Qt 
At  “  HpCp 


(D 


where  QT  is  the  total  heat  flux  (W  rrT2),  H  is 
the  thickness  of  the  water  column 
(10.7  m ),p  is  a  reference  seawater  density, 
and  Cp  is  the  specific  heat  capacity  for 
seawater,  the  increase  in  temperature  per 
unit  time  implies  a  total  heat  flux  into  the 
surface  layer  of  121 .6  WrrT2.  This  S2  net 
heat  flux  and  rate  of  temperature  increase 
is  approximately  double  that  of  SI  for  the 
same  time  period  (SI:  0.12°Cd~'). 


For  each  respective  simulation,  a  total  heat 
budget  for  the  upper  10.7  m  is  given  as 


Qt  —  Qsw(R)  +  Ql  +  Qs  +  Qnl 

+  Overt  (2) 


where  QT  is  the  total  heat  flux  (W  rrT2)  implied  by  the  change  in  the  upper  layer  temperature,  Qsw  is  the 
penetrating  solar  shortwave  flux,  R  is  the  optical  retention  of  solar  shortwave  at  10.7  m  (or  one  minus  the 
transmittance),  QL  is  the  latent  heat  flux,  Qs  is  the  sensible  heat  flux,  QNL  is  the  net  longwave  ocean-atmosphere 
exchange,  and  finally  Qvert  accounts  largely  for  the  vertical  diffusion  of  heat  out  of  the  surface  layer  but  may  also 
include  any  vertical  advective  flux.  The  convention  wherein  fluxes  out  of  the  surface  ocean  layer  are  negative 
(a  heat  loss)  is  applied.  Since  the  spatial  means  are  taken  from  a  large  horizontal  area  (41 1  km2)  in  northern  MB, 
horizontal  advection  and  diffusion  are  dropped  from  the  budget.  The  first  four  terms  on  the  right-hand  side 
(RHS)  of  equation  (2)  (the  turbulent  and  radiative  heat  fluxes  between  the  ocean  and  atmosphere;  referred  to 
herein  as  the  "direct"  fluxes)  are  recorded  in  the  modeling  system  output  at  hourly  intervals.  The  final  RFHS  term, 
the  diffusive  flux  out  of  the  surface  layer,  is  used  as  a  closure  term  for  the  heat  budget. 

Accordingly,  the  direct  fluxes  during  the  above  referenced  period  of  warming  for  S2  average  175.8  W  itT  2 
(Table  1),  and  this  implies  a  vertically  diffusive  heat  removal  (Qvert)  of  —  57  WrrT2.  In  contrast,  the  same 
analysis  for  SI  yields  a  diffusive  loss  of —75  W  rrT2  (Table  1).  Note  that  the  main  difference  in  the  direct  fluxes 
between  SI  and  S2  is  due  to  the  optical  differences  (R).  SI  relies  on  the  Paulson  and  Simpson  [1977] 
adaptation  of  the  Jerlov  IA  oceanic  water  type:  this  yields  an  R  value  of  0.77  at  1 0.7  m.  S2  relies  on  the 
biological  model  to  determine  the  attenuation  (biothermal  feedback),  and  this  R  value  at  10.7  m  averages 
0.95  during  this  period.  Based  strictly  upon  the  differences  in  penetrating  shortwave  and  allowing  for  the 


Table  1.  SI  and  S2  Fleat  Budgets:  Term-by-Term  Comparison 


S2 

SI 

(S2-S1) 

Net  heating  rate  (°C  d_1) 

0.24 

0.12 

0.12 

Net  heat  flux  (Qj;  W  m~2) 

121.6 

60.8 

60.8a 

Shortwave  attenuation  (R)  at  10.7  m 

0.96 

0.77 

0.18 

Attenuated  shortwave  (Qsw(R);  W  m~2) 

266.9 

214.1 

52.8 

Latent  heat  flux  (Ql;  W  m~2) 

-18.0 

-11.2 

-6.8 

Sensible  heat  flux  (Qs;  W  m~2) 

-1.7 

-1.8 

0.1 

Longwave  heat  flux  (Qnl;  W  m~2) 

-71.4 

-61.3 

-10.1 

Sum  of  direct  fluxes  (W  m~2) 

175.8 

139.8 

36.0b 

Implied  diffusive  flux  (Qvert;  W  m~2) 

-54.2 

-79.0 

24. 8C 

Heat  flux  below  1 0.7  md  (W  m~2) 

65.3 

142.9 

-77.6 

aNet  heat  flux  difference. 
b59.2  %  of  difference. 
c40.8  %  of  difference. 

“Convention  is  positive  for  heat  penetration. 
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increase  in  S2  latent  and  longwave 
removal  (due  principally  to  a  comparative 
increase  in  S2  SST),  the  S2  surface  layer 
should  warm  at  a  rate  exceeding  that  of  SI 
by  only  0.07°Cd~\  However,  the 
difference  in  apparent  diffusive  heat 
fluxes  out  of  the  surface  layer  increases 
this  heating  rate  difference  by  71%  up  to 
0.1 2°C  d_1;  thus,  the  S2  surface  layer  is  far 
more  retentive  than  SI. 

During  June,  S2  is  biased  warmer  than  SI 
because  periods  favoring  thermal 
stratification  in  the  upper  ocean  are  of 
longer  duration  than  periods  that  promote 
turbulent  vertical  mixing.  This  tendency 
toward  thermal  stratification  is  typical  of 
midlatitude  coastal  areas  during  the  spring 
transition  from  winter  well-mixed  water 
column  conditions  to  the  characteristic 
thermal  stratification  of  summer.  However, 
note  that  during  the  mixing  periods,  S2 
does  indeed  appear  to  lose  heat  at  a  faster 
rate  than  SI  (Figure  1 2a).  Both  the  diffusive 
flux  and  the  shortwave  transmittance 
through  10.7  m  are  greater  in  SI  than  in  S2 
(Table  1,  last  row).  This  difference 
introduces  a  cold  bias  in  the  subsurface 
waters  of  S2  (Figure  1 2b).  Hence,  when  the 
water  column  does  indeed  overturn  and 
turbulent  forces  prevail  over  the  late  spring 
and  summer  tendencies  toward  thermal 
stratification,  colder  waters  mix  into  the 
surface  layers  of  S2.  This  accounts  for  the 
reversal  in  sign  for  SST(S2  —  SI)  observed 
at  the  very  beginning  of  the  time  series  on  2  June  and  the  dramatic  lowering  of  the  SST  difference  again  on  27 
June  (Figures  5c  and  12a).  Another  consequence  of  this  apparent  biothermal  reorganization  of  thermal  energy 
fluxes  in  the  upper  ocean  is  a  simulated  increase  of  the  temporal  SST  variance:  warming  periods  promote 
comparatively  warmer  SSTs,  and  mixing  periods  promote  a  comparatively  more  rapid  SST  cooling. 

4.3.  Observations 

Are  these  simulated  differences  also  seen  in  data?  Since  we  are  unable  to  replicate  the  numerical  experiment  in 
the  "real  world"  by  removing  the  bio-optical  influence  of  phytoplankton  from  Monterey  Bay,  the  focus  of  data 
analysis  was  instead  placed  on  surface  material  property  gradients  and  comparative  intervariate  relationships. 
Accordingly,  the  chlorophyll  to  temperature  relationships  as  resolved  by  the  simulations  and  recorded  in  the 
data  were  examined.  Both  of  the  simulations  compare  favorably  to  bottle  Conductivity-Temperature-Depth 
(CTD)  data  (Figure  13)  with  respect  to  overall  changes  in  chlorophyll  magnitude  and  temperature  increases; 
however,  these  discrete  depth  interval  bottle  data  lack  the  vertical  resolution  necessary  to  determine  if  the 
differences  in  the  respective  simulations  were  also  observed.  Specifically,  the  upper  2.5  m  of  the  water 
column  in  S2  resolve  an  additional  ~0.2  units  of  normalized  temperature  increase  (Tz/T50)  that  is  not  apparent  in 
SI  (Figure  13b,  blue  shaded  portion).  Under  the  conditions  of  reduced  winds,  it  is  precisely  this  upper  ocean 
difference  in  thermal  stratification  that  is  the  main  point  of  divergence  between  the  two  simulations. 

Much  finer  vertical  scale  temperature  data  were  obtained  from  the  ScanFish  towed  platform  (SF).  This  is  a 
piloted  vehicle  that  records  data  from  an  instrument  suite  as  it  transits  from  near  surface  to  -  80  m  depth 


(A)  Mean  10.7  m  Temperature  Difference 


Figure  1 2.  (a)  Difference  (S2  —  SI )  of  the  upper  1 0.7  m  ocean  tempera¬ 
ture.  (b)  Difference  (S2  -  SI)  of  the  respective  ocean  temperature  at 
30.9  m.  All  variables  are  extracted  from  the  hourly  model  output 
increment  and  averaged  over  the  spatial  area  of  northern  MB.  The 
time  series  are  smoothed  with  a  24  h  moving  average  to  remove 
high-frequency  variability. 
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behind  the  ship.  The  data  sampling 
resolution  was  roughly  2-3  observations 
per  vertical  meter  [Jolliffet  al.,  201 2b].  An 
initial  SF  track  was  obtained  during  the 
onset  of  the  mid-June  wind  relaxation 
event  and  over  the  northwest 
continental  shelf  between  the  50  and 
200  m  isobaths  (Figure  14a).  Two  profiles 
from  the  transect  were  extracted  near 
the  50  m  isobaths  as  prototypical 
examples  of  the  temperature  and  optical 
profiles  observed  in  northern  Monterey 
Bay  during  the  field  study  (Figures  14b 
and  14c).  The  temperature  and  the  total 
absorption  coefficient  profiles  (Or(488)) 
obtained  from  the  AC-9  optical  profiling 
instrument  are  highly  correlated  with 
one  another.  Continuous  thermal 
stratification  is  evident  from  the  surface 
down  to  a  depth  of- 22-25  m.  In  both 
observed  profiles,  intense  near-surface 
(<  10  m)  thermal  gradients  on  the  order 
of  -  0.23°C  per  vertical  meter  are  evident. 
The  optical  observations  were  used  to 
calculate  the  trapping  depth,  i.e.,  the 
depth  where  -95%  of  the  total 
penetrating  solar  shortwave  irradiance 
(~  350-2500  nm)  has  been  attenuated 
[Lee  et  al.,  2005]  (see  section  A2  for 
calculation  details).  Due  to  the  very  high 
phytoplankton  biomass  in  the  surface 
layer  (independently  confirmed  by  FHigh- 
Performance  Liquid  Chromatography 
(FHPLC)  and  fluorometer  observations  during  the  field  study),  the  very  high  attenuation  of  the  photosynthetically 
available  radiation  (PAR)  spectral  band  (350-700  nm)  shoals  the  total  shortwave  trapping  depth  to 
approximately  5  m  in  both  observed  profiles.  For  comparison,  the  depths  where  95%  of  the  more  restricted  PAR 
spectral  band  has  been  attenuated  were  also  calculated:  approximately  7.5  m  for  both  observed  profiles.  The 
total  shortwave  attenuation  is  calculated  as  the  sum  of  the  infrared  (IR)  spectral  band  and  the  PAR  spectral  band. 
Attenuation  of  IR  (>  700  nm)  is  largely  invariant  due  to  the  absorption  by  pure  seawater.  PAR  attenuation,  in 
comparison,  is  severely  impacted  and  increased  by  the  presence  of  near-surface  phytoplankton. 

The  observed  temperature  data  were  used  to  estimate  the  nitrate  concentration  profiles  (*N03)  using  the 
linear  temperature-to-nitrate  relationship  observed  during  the  June  2008  MB  field  study.  Nitrate 
consumption  in  the  surface  layer  is  thus  inferred  from  the  temperature  gradient  as  well  as  the  sharp  increase 
in  aT(488).  Whereas  the  SF  did  not  measure  phytoplankton  biomass  directly,  a  separate  analysis  suggests  that 
bulk  changes  in  visible  absorption  coefficients  in  this  system  are  primarily  driven  by  changes  in  the 
concentration  of  phytoplankton  pigments  [Jolliffet  al.,  2012b]. 

Simulated  profiles  for  temperature,  nitrate,  and  phytoplankton  biomass  were  extracted  from  the 
spatiotemporally  corresponding  simulation  results  (Figures  14d  and  14e).The  sharp  near-surface  gradients  in 
observed  temperature  are  not  replicated  by  SI  (Figure  13d).  Instead,  a  much  more  gradual  temperature 
gradient  is  simulated  from  -10  to  30  m  depth.  The  simulated  phytoplankton  biomass  and  nitrate  profiles 
correspond  to  the  more  gradual  density  stratification.  It  is  important  to  note  the  inappropriate  discrepancy 
between  the  95%  PAR  attenuation  level  and  the  trapping  depth  (Figure  14d).  The  former  is  a  function  of  the 
simulated  bio-optical  attenuation,  whereas  the  latter,  serving  as  a  direct  input  to  the  physical  heating  rate 


100.00 


10.00 


mg 

Chi  1.00 
m-3 


0.10 


0.01 

1.0  1.2  1.4  1.6  1.8  2.0 


100.00 


10.00 


mg 

Chi  1.00 

m-3 


0.10 


0.01 


Figure  1 3.  A  comparison  between  simulated  and  observed  (red  asterisks) 
chlorophyll  versus  temperature  relationships.  The  temperature  is  normalized 
by  the  temperature  at  50  m  depth  (Tz/T50).  Simulations  (a)  SI  and  (b)  S2;  the 
in  situ  data  are  the  same  in  both  plots.  The  blue-shaded  portion  of  Figure  1 3b 
is  indicative  of  simulation  data  from  the  upper  2.5  m  of  the  water  column. 
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Figure  14.  (a)  Map  of  Monterey  Bay  with  ScanFish  transect  from  50  to  200  m  isobaths  indicated;  (b)  observed  profile  of  temperature,  estimated  nitrate,  and  optical 
absorption  from  the  SF  transect  near  the  50  m  isobaths.  The  dot-dashed  horizontal  line  is  the  estimated  trapping  depth;  the  thin  gray  line  is  the  depth  of  95%  PAR 
attenuation,  (c)  As  in  Figure  14b  for  the  adjacent  profile,  (d)  The  SI  simulated  profiles  of  temperature,  nitrate,  and  phytoplankton  biomass,  and  (e)  as  in  Figure  14d  for  S2. 


computations  of  the  ocean  circulation  model,  is  instead  prescribed  by  the  default  attenuation  values  of  the 
Paulson  and  Simpson  [1 977]  scheme.  These  prescribed  values  force  the  trapping  depth  in  SI  to  be  invariantly 
located  at  41  m.  This  is  grossly  in  error  when  compared  to  the  estimates  derived  from  the  optical 
observations.  As  a  result,  SI  fails  to  reproduce  the  sharp  temperature  (and  other  property)  gradients 
observed  in  the  upper  1 0  m  of  the  water  column. 

In  contrast,  the  upper  10  m  fine  structure  of  ocean  temperature  in  S2  is  closer  to  the  observations 
(Figure  14e).  The  S2  trapping  depth  is  within  8  m  of  the  observed  value  (Figures  14d).This  biologically 
enhanced  optical  attenuation  is  then  directly  communicated  to  the  physical  model  that  subsequently 
renders  an  S2  temperature  gradient  more  consistent  with  the  observations  (>  0.2oCm_1).  In  both 
simulations,  nitrate  is  above  20  mmol  N  rrT3  and  is  thus  not  considered  limiting  to  phytoplankton  growth  in 
the  biological  model  (see  section  A1,  equation  (A5)).  PAR  penetration  in  both  simulations,  in  comparison,  is 
largely  restricted  to  the  upper  20  m  of  the  water  column.  Thus,  the  additional  vertical  stability  within  the  S2 
surface  layer,  evidenced  by  the  thermal  gradient,  facilitates  an  increased  consumption  of  surface  nitrate  and  a 
corresponding  increase  in  near-surface  phytoplankton  biomass. 
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These  simulated  profile  differences  underscore  the  model  result  contrasts  analyzed  in  previous  sections. 
However,  herein  it  is  shown  that  S2  is  more  successful  at  reproducing  the  observed  near-surface  gradients 
and  trends  in  temperature,  nitrate,  and  potentially  phytoplankton  biomass,  insofar  as  the  optical 
observations  may  serve  as  a  biomass  proxy.  The  combined  observations  and  simulations  suggest  that  light 
availability,  rather  than  macronutrient  concentration,  is  the  main  constraint  upon  phytoplankton  growth 
within  this  regime.  Hence,  the  vertically  stabilizing  effect  of  biothermal  feedback  tends  to  concentrate 
phytoplankton  biomass  near  the  surface  where  nutrients  may  be  utilized  for  growth  more  efficiently.  This 
mechanism,  in  turn,  shoals  the  trapping  depth  nearer  the  surface. 

5.  Discussion  and  Conclusions 

Ramp  et  at.  [1991]  observed  that  under  conditions  of  low  wind  stress  and  high  solar  insolation  in  coastal 
California,  large  vertical  gradients  in  upper  ocean  temperatures  occurred  in  patches  consistent  with  high  values 
of  surface  chlorophyll.  They  suggest  that  on  the  basis  of  their  observations  and  simple  one-dimensional  heating 
rate  calculations,  it  is  the  additional  optical  attenuation  provided  by  surface  phytoplankton  stocks  that  brings 
about  the  observed  surface  heating  and  associated  temperature  gradient.  The  results  presented  here  and  in 
earlier  work  [Jolliffet  a!.,  2012a]  support  this  hypothesis.  However,  here  it  is  shown  that  the  additional  optical 
attenuation  may  not  simply  impact  the  vertical  temperature  structure  and  then  cease  to  have  any  further 
impact;  rather,  the  vertical  temperature  structure  then  feeds  back  to  the  biology  to  concentrate  biomass  near 
the  surface  and  improve  the  nutrient  utilization  for  growth.  To  use  an  engineering  frame  of  reference,  this  is  an 
example  of  a  positive  feedback  loop. 

The  observations  of  Ramp  et  al.  [1991]  indicate  that  this  biophysical  feedback  loop  may  reach  extremes  we 
are  unable  to  presently  simulate  with  our  model:  a  4.7°C  gradient  in  the  upper  2  m  of  ocean.  Such  fine-scale 
gradients  are  below  the  vertical  grid  scale  of  many  ocean  models.  Nevertheless,  our  two-way  coupled  ocean- 
atmosphere  simulations  strongly  suggest  that  such  situations  would  have  a  profound  impact  on  air-sea 
exchange  of  thermal  energy  and  potentially  many  other  upper  ocean  biological  and  physical  processes.  For 
example,  an  enhanced  near-surface  stratification  and  associated  reduction  in  vertical  diffusion  may  have 
important  consequences  for  air-sea  gas  exchange  that  may  impact  carbon  and  oxygen  budgets. 

Vertical  diffusion  of  ocean  scalar  quantities  in  this  modeling  system  is  quantified  via  the  Mellor  Yamada  2.5 
vertical  turbulence  closure  scheme  [Mellor  and  Yamada,  1982].  This  numerical  recipe  is  one  approximation 
out  of  many  other  possible  choices  [e.g.,  Large  et  al.,  1 994],  Any  numerical  approximation  for  vertical  mixing 
of  scalars  is  likely  to  be  imperfect.  Nevertheless,  the  observations  of  Ramp  etal.  [1991]  support  the  hypothesis 
that  under  the  conditions  of  low  wind  stress  and  high  solar  insolation,  surface  phytoplankton  blooms  provide 
a  biothermal  feedback  effect  that  may  suppress  vertical  diffusion  of  thermal  energy  and  other  material 
properties.  In  the  particular  circumstance  of  Monterey  Bay,  California  during  the  upwelling  season,  such 
increases  in  vertical  stability  are  beneficial  to  the  flora  due  to  the  relative  abundance  of  macronutrients  and 
the  overall  instability  of  this  highly  turbulent  hydrodynamic  regime. 

The  concept  of  flora  having  an  impact  on  the  physical  environment  that,  in  turn,  further  promotes  the  growth 
of  the  flora  is  analogous  to  the  concept  of  "microclimates"  in  terrestrial  vegetation.  The  most  obvious  example 
is  the  shading  effect  of  a  forest  canopy:  reducing  incident  PAR  to  sublethal  levels,  increasing  moisture  retention, 
and  reducing  temperatures  and  respiration  costs  [Callaway  and  Pugnaire,  2007,  and  references  therein].  All 
of  these  effects  are  a  perceived  benefit  to  many  species  of  vegetation  emerging  on  the  forest  floor.  Discovering 
a  marine  analogue  of  this  situation  for  phytoplankton  may  at  first  seem  dubious.  The  phytoplankton  are,  to 
first  order,  subject  to  the  comparatively  extreme  physical  disturbance  and  energy  fluctuations  of  the  surface 
ocean  and  are  little  more  than  passive  recipients  of  whatever  physical  conditions  may  exist  there.  However,  the 
evidence  presented  here  and  elsewhere  shows  that  under  some  specific  local  conditions  of  reduced  wind 
stress,  high  solar  insolation,  and  available  macronutrients,  the  optical  properties  of  phytoplankton  may  indeed 
create  a  beneficial  "microclimate"  of  vertical  stability  in  the  upper  ocean  (<  10  m  depth)  that  promotes  the 
growth  and  persistence  of  a  surface  phytoplankton  bloom  (>  2  mg  rrT3  surface  chlorophyll). 

This  biothermal  potential  of  phytoplankton  must  also  be  considered  with  respect  to  diurnal  warming.  The 
diurnal  thermocline  and  features  of  vertical  stratification  may  be  exploited  by  phytoplankton  if  sufficient 
nutrients  are  available.  In  a  separate  study ,  Jolliffet  al.  [2012b]  examined  the  apparent  fluorescence  efficiency 
of  phytoplankton  from  collocated  fluorescence  and  optical  data.  The  persistent  surface  warm  layer  near  the 
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coasts  within  Monterey  Bay  was  certainly  exploited  by  diatoms  (identified  there  by  HPLC  data)  growing  at 
eutrophic  concentrations  (>2  mg  chi  rrT3).  Away  from  the  coasts,  however,  the  apparent  fluorescence 
efficiency  within  an  obviously  diurnal  surface  warm  layer  suggested  an  impairment  of  the  photosynthetic 
apparatus  near  the  surface  and  chlorophyll  observations  there  were  generally  below  -  2  mg  rrT3.  It  was 
hypothesized  that  the  availability  of  dissolved  iron  may  impact  the  ability  of  phytoplankton  in  general,  and 
diatoms  in  particular,  to  adapt  to  high-light  conditions  or  highly  variable  light  conditions  near  the  air-sea 
interface  and  thus  exploit  the  diurnal  warm  layer.  The  biological  model  integrated  into  COAMPS  herein  did  not 
consider  dissolved  iron  as  a  limiting  nutrient  or  explicit  state  variable.  Future  work  may  have  to  consider  iron 
dynamics  more  explicitly,  or  other  biogeochemical  processes  that  impact  phytoplankton  growth  efficiency. 

In  conclusion,  as  integrated  model  simulations  of  the  ocean-atmosphere  system  move  toward  finer  scales  of 
spatial  and  temporal  resolution,  investigators  will  need  to  consider  the  role  of  the  biosphere  in  the  ocean- 
atmosphere  exchanges  of  energy  and  all  of  the  processes  subsequently  impacted  by  this  exchange.  More 
broadly,  the  practice  of  modeling  ocean  physics  independently  of  ocean  biology  and  vice  versa  may  have  to 
cease.  The  situation  is  similar  to  that  of  meteorology  and  physical  oceanography:  one  system  (the  ocean) 
cannot  simply  be  prescribed  as  an  invariant  boundary  condition  to  the  other  (the  atmosphere).  Instead,  the 
dynamical  interactions  between  the  atmosphere  and  hydrosphere  must  be  accounted  for  within  the 
simulation  itself.  Our  results  strongly  suggest  that  simulations  of  the  upper  ocean  temperature  structure, 
ocean-atmosphere  exchange  of  thermal  energy,  and  some  aspects  of  the  wind  field  and  surface  ocean 
circulation  are  dynamically  sensitive  to  biothermal  effects,  and  so  the  potential  dynamical  influence  of  the 
biosphere  must  also  be  considered. 

Appendix  A:  Model  Formulae 
A1.  Biological  Model 

The  biological  model  integrated  into  COAMPS  is  a  four-component  nutrient-phytoplankton-detritus  (NPD) 
model.  Nitrate  is  initialized  from  observed  temperature-to-nitrate  relationships  [Chavez  et  al.,  1996],  and 
subsequent  nitrogen  cycling  is  allowed  to  "spin-up"  as  the  simulated  phytoplankton  biomass  assimilates 
nitrogen  and  is  then  subject  to  implicit  grazing/mortality  terms.  Nitrogen-containing  particulate  detritus  sinks 
and  is  subject  to  respiration.  Mortality  and  respiration  contribute  to  the  ammonium  pool,  and  the  cycle  is 
completed  via  light-inhibited  nitrification.  There  is  presently  no  benthic  component  in  the  model;  particulate 
detritus  that  reaches  the  bottommost  grid  cell  is  respired  back  to  the  ammonium  pool.  Conversion  between 
phytoplankton  carbon  and  nitrogen  is  accomplished  via  the  Redfield  ratio  [Redfield  et  al.,  1963]. 

The  conservation  equation  for  phytoplankton  carbon  (C,  mmol  C  m~3)  follows  the  NCOM  conservation 
equation  for  temperature  given  in  Barron  et  al.  [2006]: 

^=-V-(\ /C)+A„V3C  +  !(k„J0+F(C)  (Al) 

Advection  (first  term  on  the  right-hand  side  (RFHS)),  horizontal  diffusion  (second  term  on  RHS),  and  vertical 
diffusion  (third  term  on  RHS)  are  solved  via  NCOM.  The  last  term  on  the  RHS  represents  the  biological  sources 
and  sinks  that  are  described  in  detail  herein.  Equation  (Al)  is  applicable  to  the  other  three  biological  model 
state  variables  that  are  quantified  with  respect  to  nitrogen:  nitrate  (A/,,  mmol  N  itT3),  ammonium  (A/a, 
mmol  N  itT3),  and  particulate  detritus  (D,  mmol  N  itT3).  A  complete  list  of  parameters,  state  variables,  and 
internal  variables  for  the  biological  model  is  given  in  Table  Al .  The  internal  biological  model  time  step  (At)  is 
5  s.  Each  state  variable  is  defined  over  the  entire  three-dimensional  coordinate  domain  ( x ,  y,  z).  Below,  these 
spatial  subscripts  are  dropped  for  simplicity,  except  where  required  to  indicate  depth  (z). 

Phytoplankton  growth  is  a  function  of  available  light  and  nitrogen.  Light  is  quantified  as  downwelling 
photosynthetically  available  radiation  ( EPAr ,  W  rrT2).  Description  of  the  light  attenuation  is  provided  in 
section  A2  as  part  of  the  optical  computations.  Expanding  here  the  last  term  on  the  REIS  of  equation  (Al ): 

m  =^  =  (Mr-  rr)C  (A2) 

The  realized  rate  of  growth  (//r)  is  the  minimum  of  the  light-limited  (,kl)  and  nitrogen-limited  (//N)  growth 
rate  calculations: 
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Table  A1.  Biological  Model  Parameters,  State  Variables,  and  Internal  Variables 


Parameter 

Description 

Units 

Value 

Qn 

Nitrogen  half-saturation  constant 

mmol  N  rrF3 

1.0 

A  max 

Maximum  phytoplankton  growth  rate 

s_1 

3.2  X  1 0  ! 

Iv 

Ivlev  grazing  parameter 

(mmol  C  rrF3)-1 

5.56 

^max 

Maximum  phytoplankton  grazing  rate 

s-1 

5.8  X  1 0' 

Si 

PAR  saturation  constant 

WnT2 

145 

Basal  mortality  rate 

s-1 

1 .2  X  1 01 

¥>l 

Grazing  fraction-to-detritus  pool 

— 

0.50 

* 

Molar  nitrogen-to-carbon  ratio 

mol  moF1 

0.1509 

ws 

Detritus  sinking  rate 

m  s-1 

2.3  X  1 0 

Bl 

Detritus  loss  rate 

s1 

2.3  X  1 0' 

Qni 

Nitrification  half-saturation  constant 

mmol  N  rrF3 

1.0 

r2 

Maximum  nitrification  rate 

s-1 

8.2  X  1 0  ! 

Ct 

C  threshold  concentration 

mmol  C  m~3 

1 .0  X  1 0~' 

State  variable 

C 

Phytoplankton  carbon 

mmol  C  m~ 

N\ 

Nitrate  (-mitrite) 

mmol  N  itF3 

Na 

Ammonium 

mmol  N  rrF3 

D 

Particulate  detritus 

mmol  N  rrF3 

Internal  variable 

Ar 

Realized  phytoplankton  growth  rate 

s-1 

AL 

Light-limited  phytoplankton  growth  rate 

s-1 

An 

Nitrogen-limited  phytoplankton  growth  rate 

s-1 

Tr 

Realized  phytoplankton  loss  rate 

s-1 

A/t 

dissolved  inorganic  nitrogen 

mmol  N  rrF3 

NX 

/VA  removal  flux  via  nitrification 

mmol  N  rrF3  s~' 

Ar  =  MINCul^n) 


(A3) 


where  the  light-limited  growth  rate  is  a  function  of  £PAR  calculated  for  vertical  level  (z),  the  PAR  saturation 
constant  (Si),  and  the  maximum  growth  rate  (/tmax,  Table  A1)  [cf.  Webb  et  a!.,  1974]: 


Al=Amax(i  -e-[£RA»(*n) 


Nitrogen  limitation  is  a  function  of  the  total  nitrogen  available  (A/t  =  A/,  +  A/a): 

A/j 


An  ~  Amax 


(A/t  +  On) 


(A4) 


(A5) 


The  realized  rate  of  grazing  is  based  on  the  Ivlev  formulation  [Ivlev,  1955;  Hallam,  1978;  Franks  eta!.,  1986]  and 
requires  a  maximum  potential  rate  (Tm)  and  an  Ivlev  parameter  (Iv).  The  herbivore  grazing  rate  saturates 
toward  the  maximum  rate  at  elevated  phytoplankton  biomass  (C).  Here  a  basal  mortality  rate  (Bm)  is  added  to 
the  expression: 

-lvC\ 


rr  =  rm  (i  -e 


firr 


(A6) 


If  the  phytoplankton  carbon  falls  below  a  threshold  value  (Ct)  then  Tr  is  set  to  0. 

Total  grazing  and  mortality  of  phytoplankton  biomass  provides  a  source  to  the  remaining  nitrogen 
compartments.  The  biological  terms  for  particulate  detritus  (D)  are  given  as  follows: 

AD  ,  /AD\ 

F(D)  =  ^7  =  rr  C/f'  ~R'D 


(A7) 


where  the  grazed  phytoplankton  carbon  is  converted  to  nitrogen  (x)  and  multiplied  by  the  fraction  of  grazed 
nitrogen  allotted  to  the  D  pool  Sinking  is  calculated  from  the  sinking  speed  parameter  (ws)  and  the 
vertical  detritus  gradient.  The  loss  rate  (/?i)  represents  the  implicit  bacterial  respiration  back  to  the 
ammonium  compartment  (/VA).  If  the  depth  (z)  is  at  the  maximum  depth  increment  then  D(z)  is  returned  to 
the  Na(z)  pool. 
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Table  A2. 

Symbol 

Optical  Model  Parameters 

Description 

Units 

Value 

6 

Carbon-to-chlorophyll  ratio  (mass) 

gg-1 

47.2 

0C-] 

Empirical  absorption  model  coefficient  1 

— 

-0.017 

a2 

Empirical  absorption  model  coefficient  2 

— 

0.076 

a3 

Empirical  absorption  model  coefficient  3 

— 

0.554 

«4 

Empirical  absorption  model  coefficient  4 

— 

-1.149 

Ai 

Empirical  scattering  model  coefficient  1 

— 

0.416 

h 

Empirical  scattering  model  coefficient  2 

— 

0.766 

bbw(490) 

Pure  seawater  backscattering  coefficient 

m-1 

0.0012 

f2 

Fraction  of  shortwave  in  PAR  spectrum 

— 

0.48 

£l 

K-\  PAR  attenuation  coefficient  1 

m-1 

-0.057 

£ 2 

K-\  PAR  attenuation  coefficient  2 

— 

0.482 

£ 3 

Kf  PAR  attenuation  coefficient  3 

— 

4.221 

K2  PAR  attenuation  coefficient  1 

m-1 

0.183 

52 

l<2  PAR  attenuation  coefficient  2 

— 

0.702 

?3 

K2  PAR  attenuation  coefficient  3 

— 

-2.567 

Accordingly,  the  ammonium  (A/a)  biological  terms  are  given  as  follows: 

~^f  =  Tr  CxO-fi)+RiD-p,CXj±-NX 

where  the  last  term  on  the  RHS  (A/t)  is  the  ammonium  mass  flux  lost  to  nitrification: 


/Va  =  A/a 


A/a 


Fb 


V,A/A  +  Qn 

Nitrification  is  only  simulated  in  darkness:  if  EPAR(z)  >  0,  then  A/a(z)  =  0. 


(A8) 


(A9) 


The  model  does  not  distinguish  nitrate  versus  nitrite  oxidation  states;  thus,  the  nitrate  state  variable's 
biological  terms  (A/|  =  nitrate  +  nitrite)  are  given  as  follows: 


A  A/| 


=  A/a  —  MrQ 


A/| 


At  H  ‘  ~  Nt 

A/,  is  initialized  from  empirical  temperature-to-nitrate  relationships,  and  this  is  maintained  at  the  model 
domain's  boundary  value  points. 


(A10) 


At  4  components  and  13  parameters,  this  biological  model  is  intended  to  present  a  minimal  level  of 
ecosystem  and  biogeochemical  complexity.  Phytoplankton  functional  groups,  explicit  zooplankton  biomass, 
other  nutrient  sources/requirements,  dissolved  nonliving  organic  matter,  a  benthic  component,  and  an 
explicit  microbial  loop  are  all  omitted  from  this  initial  formulation.  This  model  provides  a  framework  wherein 
these  additional  terms  and  sets  of  equations  may  be  added,  as  required,  to  improve  model  performance. 


A2.  Optical  Computations 


Biological  feedback  to  the  physics  is  represented  in  the  system  via  the  variable  penetration  of  solar  shortwave 
energy  into  the  surface  ocean  that  is  based  upon  the  optical  attenuation  estimated  in  the  biological  model.  This 
interface  between  the  biological  and  physical  ocean  models  replaces  the  "Jerlov"  water  types  [Jerlov,  1976] 
previously  used  in  the  NCOM  heating  rate  computations  [Jolliffet  al.,  201 2a].  The  attenuation  of  solar  shortwave 
in  the  PAR  spectral  range  (350-700  nm)  is  calculated  identically  in  the  physical  and  biological  model 
components.  This  calculation  is  based  upon  Lee  etal.  [2005].  This  method  requires  an  estimate  of  the  surface  total 
absorption  coefficient  at  490  nm  [aT  (490)]  and  the  total  surface  backscattering  coefficient  at  490  nm  [bbT  (490)]. 
These  quantities  are  calculated  from  the  surface  phytoplankton  carbon  (Q  biological  state  variable  as  follows: 

logi0[aT(490)]  =  a-|dog3  +  a2clog2  +  a3clog  +  a4  (All) 


and 

clog  =  log10  (C*6*  1 )  (A12) 

where  C*  is  the  phytoplankton  carbon  in  mass  units  (mg  C  rrT3)  and  6  is  the  carbon-to-chlorophyll  mass  ratio 
(Table  A2).  All  optical  coefficient  values  are  given  in  Table  A2.  Equation  (All)  is  an  empirical  relationship 
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between  phytoplankton  chlorophyll  and  total  absorption  coefficients  obtained  from  in  situ  data  from 
Monterey  Bay  [ Jolliffet  al.,  2012b].  The  carbon-to-chlorophyll  ratio  value  (Table  A2)  is  also  obtained  from 
these  data.  This  empirical  expression  implicitly  includes  the  absorption  signal  that  may  arise  from  nonliving 
organic  matter  associated  with  phytoplankton  productivity. 

The  calculation  of  the  total  backscattering  coefficient  [bbT  (490)]  is  based  on  literature.  The  particle  scattering 
coefficient  at  550  nm  [bp(550)]  as  a  function  of  chlorophyll  concentration  is  taken  from  Morel  and 
Maritorena  [2001]: 

[bp(550)]  =^(Cd-'Y2  (A13) 


The  total  backscattering  coefficient  at  490  nm  is  calculated  as  a  function  of  bp  (550): 

[frbT(490)]  =  jo.002  +  0.01  [0.50  -  0.25(clog)[  j  [bp(550)]  +  bbw(490) 


where 


and 


v  =  (0.5)(clog  -  0.3),  0.02  <  [C*0_1]  <  2  mg  m~3 
v  =  0,  [CfT1]  >  2  mg  rrT3 


(A14) 


(A15) 

(A16) 


Once  these  optical  quantities  are  determined  [aT  (490),  bbT(490)],  the  Lee  et  al.  [2005]  algorithm  is  employed 
to  attenuate  simulated  shortwave  into  the  surface  ocean.  Total  shortwave  (fT,  W  rrT2)  just  below  the  air-sea 
interface  (z  =  CT)  is  partitioned  into  PAR  (fPAR,  W  rrT2)  and  infrared  (f|R_  W  rrT2)  spectral  components.  This 
follows  the  Lee  et  al.  [2005]  spectral  partition  for  PAR  (350-700  nm)  and  IR  (700-2500  nm). 

EpAR(0-)  =  EToT(0-)f2,  Eir(0-)  =  £Tot((T)(1  -  f2)  (A17) 

It  is  assumed  that  48%  of  the  total  shortwave  is  in  the  PAR  spectral  range  (Table  A2).  The  attenuation 
coefficients  for  PAR  (/CPAr(z))  and  IR  (K)R(z))  are  both  a  function  of  depth: 

fPAR(z)  =  £PAR(0-)e-*PAR(z)z  (Al  8) 

and 

£ir(z)  =  £IR(0-)e-K»^  (A  1 9) 


The  depth-dependent  attenuation  coefficient  for  IR  (K|R(z))  is 


K)r(z)  =  0.56  + 


2.304 

(0.001  +z)0'65 


(A20) 


PAR  attenuation  is  also  from  Lee  et  al.  [2005]  but  requires  several  additional  computation  steps: 

/<PAR  (z)  =  K,  +  —  Q+2  0.5  (A2 1 ) 

where 

K,  =  £i  +  e2[aT(490)]°-5  +  e3  [fc>bT(490)[  (A22) 

and 

Ki  =  Ci  +  CilM 490)]  +  C3[bbT(490)]  (A23) 

The  attenuation  of  £,R  and  £PAr  are  required  for  the  NCOM  heating  rate  computations.  The  incident  total 
shortwave  is  simulated  via  the  COAMPS  atmospheric  model  component.  The  attenuation  of  £PAr  is  also 
required  for  the  biological  model:  the  value  of  PAR  at  each  vertical  level  (z)  is  used  in  equation  (A4)  to 
determine  the  light-limited  rate  of  growth. 

It  is  important  to  note  that  given  the  above  equations,  the  phytoplankton  biomass  state  variable  is  the  driver 
of  changes  in  the  total  absorption  and  total  backscattering  coefficients,  and  the  subsequent  alterations  to 
/<pAr(z).  This  Case  I -type  assumption  [Morel  and  Prieur,  1977]  is  appropriate  for  Monterey  Bay,  California, 
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and  other  marine  environments  wherein  phytoplankton  and  associated  nonliving  organic  matter  are  the 
dominant  bio-optical  constituents.  Two  important  points  need  top  be  clarified  with  respect  to  the 
biological  and  optical  formulations  presented  herein.  First,  we  recognize  and  caution  that  other  processes 
such  as  terrestrial  runoff,  dispersion  of  river  plumes,  and  suspension  of  mineral  sediments  (i.e..  Case  2) 
would  violate  our  optical  assumptions.  These  additional  processes  will  be  addressed  in  future  work  by  the 
addition  of  terms  to  the  state  variable  space  and  associated  optical  terms  to  equations  (A1 1 )  and  (A1 4). 

Second,  the  sinking  mass  of  organic  detritus  represented  explicitly  in  our  biological  state  variable  space 
(D,  equation  (A7))  is  not  represented  by  a  specific  mass  to  optical  property  mathematical  transform  in  our 
optical  equations.  Instead,  equations  (A1 1)  and  (A14)  use  the  Case  1  assumption  to  render  the  total 
absorption  and  backscattering  coefficients  as  an  empirical  function  of  the  phytoplankton  biomass.  In  the 
case  of  absorption  coefficients,  this  is  because  the  main  contributor  from  nonliving  organic  matter  is  due  to 
colored  detrital  matter,  or  CDM  [see  Siegel  et  at.,  2005],  which  includes  a  largely  dissolved  component.  These 
materials  are  generally  refractory,  depleted  in  nitrogen  [Kahlerand Koeve,  2001],  and  require  additional  sets  of 
state  variables  and  parameterizations  than  those  presented  within  our  simple  NPD  model  framework.  Future 
work  will  move  toward  more  explicit  representations  of  CDM  in  both  the  biological  and  optical  model 
equations. 
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